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Abstract 

We propose a novel theory of dark matter (DM) superfluidity that matches the successes of the ACDM 
model on cosmological scales while simultaneously reproducing the Modified Newtonian Dynamics (MOND) 
phenomenology on galactic scales. The DM and MOND components have a common origin, representing 
different phases of a single underlying substance. DM consists of axion-like particles with mass of order 
eV and strong self-interactions. The condensate has a polytropic equation of state P ~ p 3 giving rise to a 
superfluid core within galaxies. Instead of behaving as individual collisionless particles, the DM superfluid 
is more aptly described as collective excitations. Superfluid phonons, in particular, are assumed to be 
governed by a MOND-like effective action and mediate a MONDian acceleration between baryonic matter 
particles. Our framework naturally distinguishes between galaxies (where MOND is successful) and galaxy 
clusters (where MOND is not): due to the higher velocity dispersion in clusters, and correspondingly higher 
temperature, the DM in clusters is either in a mixture of superfluid and normal phase, or fully in the normal 
phase. The rich and well-studied physics of superfluidity leads to a number of observational signatures: 
array of low-density vortices in galaxies, merger dynamics that depend on the infall velocity vs phonon sound 
speed; distinct mass peaks in bullet-like cluster mergers, corresponding to superfluid and normal components; 
interference patters in super-critical mergers. Remarkably, the superfluid phonon effective theory is strikingly 
similar to that of the unitary Fermi gas, which has attracted much excitement in the cold atom community 
in recent years. The critical temperature for DM superfluidity is of order mK, comparable to known cold 
atom Bose-Einstein condensates. Identifying a precise cold atom analogue would give important insights on 
the nricrophysical interactions underlying DM superfluidity. Tantalizingly, it might open the possibility of 
simulating the properties and dynamics of galaxies in laboratory experiments. 


1 Introduction 

The most clear-cut evidence for dark matter (DM) comes from observations on the largest scales. 
The standard A-Cold-Dark-Matter (ACDM) model, in which DM consists of collisionless particles, 
does exquisitely well at fitting the background expansion history, the detailed shape of microwave 
background and matter power spectra, as well as the abundance and mass function of galaxy 
clusters. On smaller scales, however, the situation is murkier. As simulations and observations of 
galaxies have improved, a number of challenges have emerged for the CDM paradigm. 

For starters, galaxies in our universe are observed to be remarkably regular, a fact embodied by 
various empirical scaling relations. The most striking example is the Baryonic Tully Fisher Relation 
(BTFR) |TH3] , which relates the baryonic mass Mb to the asymptotic circular velocity u c 0 

M b ~u c 4 . (1) 

1 The BTFR extends the old Tully-Fisher relation [5], relating the optical luminosity to velocity as L ~ vt- 
Since the mass-to-light ratio is not constant among different types of galaxies, the inferred slope and scatter end up 
depending on the choice of band filter. Replacing luminosity by total baryonic mass PE| (i.e., stars and gas) reduces 
the scatter and extends the validity the scaling relation over many decades in mass [4], as shown in Fig. fll 
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Figure 1: The Baryonic-Tully-Fisher-Relation (BTFR), reproduced from [6]. The dark blue points are star- 
dominated galaxies; the light blue circles are gas-dominated. The dashed line has a slope of 3, corresponding 
to the ACDM prediction. The dotted line has slope 4, in good agreement with the data. 

Figure [I] reproduced from [6], shows excellent agreement with remarkably little scatter in the 
high-mass end comprised of star-dominated (dark blue circles) and gas-dominated disc galaxies 
(light blue circles). On the theory side, the standard collapse model predicts a scaling between the 
total mass (dark plus baryonic) and circular velocity at the virial radius: M v j r ~ v^ ir . Despite the 
different slope, this is not a priori inconsistent with 0 since the translation from virial parameters 
to observables can be mass-dependent. However, the real challenge for ACDM lies in explaining the 
remarkably small level of scatter around this slope in the high-nrass end, as shown in Fig. [lj How 
can baryonic feedback processes, which are inherently stochastic, result in such a tight correlation 
across different galaxy types? Indeed, recent hydrodynamical simulations [7] show considerably 
larger scatter than observations [lj. 

Another set of challenges comes from dwarf satellite galaxies in the Local Group. Dwarf satellites 
are highly DM-dominated objects and thus well-suited to detailed tests of DM microphysics. As 
the old “missing satellite” problem [Hf-fTUlj has gradually been alleviated through the discovery of 
ultra-faint dwarfs EHH, new sharper problems have emerged. Recent attempts at matching the 
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populations of simulated subhaloes and observed MW dwarf galaxies have revealed a “too big 
to fail” problem ca ns]: the most massive dark halos are too dense to host the brightest MW 
satellites. Even more puzzling is the fact that the majority of the MW P2H2Q! and Andromeda 
(M31) [211423 j satellites lie within vast planar structures and are co-rotating within these planesj^] 
This is puzzling for ACDM, though mechanisms have been proposed I25H28] . 

Another puzzle comes from tidal dwarfs — “recycled” galaxies that form in the tidal material 
created by merging spirals. While standard theory tells us that tidal dwarfs should be devoid of 
dark matter |29t431| . recent observations of three such objects around NGC5291 [32] have revealed 
a dynamical mass discrepancy of about 2-3 times the visible mass. A standard explanation is that 
the missing matter is in the form of cold baryonic gas [32] , however this seems unlikely given their 
flat rotation curves and the remarkable consistency with the BTFR [33]. One should be wary of 
drawing definitive conclusions from a few objects, but a larger sample of tidal dwarfs will reduce 
uncertainties and can provide a critical test of ACDM. 


1.1 MOND: successes, challenges and failures 


A radical alternative is Modified Newtonian Dynamics (MOND) [34H37] . which proposes to replace 
DM with a modification of the Newtonian force law. The force law is standard at large acceleration 
(a ~ on for on 3> «o) but modified at low acceleration (a ~ y / ONOo for on <C ao). This empirical 
force law has been remarkably successful at explaining a wide range of galactic phenomena [5, 53]. 
For spiral galaxies, it predicts asymptotically flat rotation curves and provides an excellent fit to 
detailed rotation curves [38] . The critical acceleration ao is the only free parameter (apart from the 
0(1) mass-to-light ratio for each galaxy), with the best-fit value intriguingly of order the present 
Hubble parameter: 

ao — -Ho ~ 1.2 x 10 -8 cm/s 2 . (2) 
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The BTFR is an exact consequence of this force law — deep in the MOND regime (on ao), a 

test particle will orbit an isolated spherically-symmetric source according to f = yJ G *% ba ° , and 
hence 


M b 


Gnoo 


(3) 


The vast planar structures seen around the MW and Andromeda also find a plausible explanation 
in MOND, as the result of tidal stripping during a fly-by encounter between these galaxies. With the 
MOND force law, this encounter has been estimated to have occurred ~ 10 Gyr ago, with < 55 kpc 
closest approach distance (39]. Unlike in ACDM, where galaxies are surrounded by extended DM 
halos and dynamical friction would cause a rapid merger, in MOND there is only stellar dynamical 
friction and a merger can be avoided [401442] . MOND predicts that tidal dwarf galaxies should have 
flat rotation curves and fall on the BTFR, consistent with the NGC5291 dwarfs [32j 33] . 

On the flip side, dwarf satellites, particularly the MW dwarf spheroidals, have long posed a 
challenge for MOND [431447] . Five of the classical dwarfs are consistent with the BTFR, but two 
(Draco and Ursa Minor) fall below it [451 144] . Nearly all the ultra-faint dwarfs lie systematically 
below the BTFR [43]. However, the derivation of the BTFR in MOND assumes dynamical equilib¬ 
rium, whereas the discrepant dwarfs may be undergoing tidal disruption [48] , Moreover, velocity 
estimates for these objects are complicated by interlopers [49]. On the other hand, MOND does an 

“Phase-space correlated dwarfs have also been found around galaxies beyond the Local Group l24l . 
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excellent job at explaining the observed velocity dispersions in Andromeda’s dwarf satellites 150. 51. 
Finally, globular clusters also pose a challenge for MOND [52] . 

MOND faces much more severe challenges on extra-galactic scales. To reproduce the observed 
temperature profile of galaxy clusters [53], one must invoke some form of dark matter, either 
as massive neutrinos [541456] and/or cold dense gas clouds [57]. Relativistic versions of MOND, 
such as the Tensor-Vector-Scalar (TeVeS) theory [58)463] and other related proposals [65)468] (see 
[69 ] for a review), cannot match the CMB power spectrum Earn]. Without a significant dark 
matter component, the baryonic oscillations in the matter power spectrum tend to be far too 
pronounced mm- Finally, numerical simulations of MONDian gravity with massive neutrinos 
fail to reproduce the observed cluster mass function mm- 

1.2 DM-MOND hybrids 

What we have learned is that MOND and CDM are each successful in almost mutually exclusive 
regimes. The ACDM model successfully explains the expansion and linear growth histories, as well 
as the abundance of clusters, but faces a number of challenges on galactic scales. MOND does very 
well overall at explaining the observed properties of galaxies, in particular the empirical scaling 
relations, but it seems highly improbable that it can ever be made consistent with the detailed 
shape of the CMB and matter power spectra. 

This has led various people to propose hybrid models that include both DM and MOND phe¬ 
nomena [751482] . For instance, one of us recently proposed such a hybrid model, involving two 
scalar fields [83]: one scalar field acts as DM, the other mediates a MOND-like force law. This 
model enjoys a number of advantages compared to TeVeS and other relativistic MOND theories. 
For starters, it only requires two scalar fields, as opposed to the scalar and vector fields of TeVeS. 
Secondly, unlike TeVeS, its predictions on cosmological scales are consistent with observations, 
thanks to the DM scalar field. Finally, the model offers a better fit to the temperature profile of 
galaxy clusters. 

The improved consistency with data does come at the price of having two a priori distinct 
components — a DM-like component and a modified-gravity component. It would be much more 
compelling if these two components somehow had a common origin. Furthermore, the theory must 
be adjusted such as to avoid co-existence of DM-like and MOND-like behavior. This requires 
that the parameters of the theory be mildly scale or mass dependent, which adds another layer of 
complexity. 

1.3 Unified approach: MOND phenomenon from DM superfluidity 

In this paper, along with its shorter companion [83], we propose a unified framework for the DM 
and MOND phenomena. The DM and MOND components have a common origin, representing 
different phases of a single underlying substance. This is achieved through the rich and well-studied 
physics of superfluidity. 

There are two central ideas underlying this work. The first idea is quite general, namely that DM 
forms a superfluid inside galaxies with a coherence length of order the size of galaxies. As we will 
see, the phenomenon of DM superfluidity is quite generic if the DM particle is sufficiently light and 
has sufficiently strong self-interaction. Specifically, as a back-of-the-envelope calculation, we can 
estimate the condition for the onset of superfluidity by ignoring interactions among DM particles. 
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With this simplifying approximation, the requirement for superfluidity amounts to demanding that 
the de Broglie wavelength AdB ~ J,, °f DM particles should overlap. Using the typical velocity v 
and density of DM particles in galaxies, this translates into an upper bound m < 2 eV on the DM 
particle mass. 

Another requirement for Bose-Einstein condensate is that DM thermalize within galaxies. We 
assume that DM particles interact through contact repulsive interactions. Demanding that the 
interaction rate be larger than the galactic dynamical time places a lower bound of 0.1 cm 2 /g. 
This is just below the most recent constraint < 0.5 cm 2 /g from galaxy cluster mergers [SS], though 
we will argue such constraints must be carefully reanalyzed in the superfluid context. 

Again ignoring interactions, the critical temperature for DM superfluidity is T c ~ mK, which 
intriguingly is comparable to known critical temperatures for cold atom gases, e.g., 'Li atoms have 
T c ~ 0.2 mK. We will see that cold atoms provide more than just a useful analogy — in many 
ways, our DM component behaves exactly like cold atoms. In cold atom experiments, atoms are 
trapped using magnetic fields; in our case, it is gravity that attracts DM particles in galaxies. 

The superfluid nature of DM dramatically changes its macroscopic behavior in galaxies. Instead 
of behaving as individual collisionless particles, the DM is more aptly described as collective exci¬ 
tations, which at low energy are just phonons. In the non-relativistic regime and at lowest order 
in derivatives, it is well-known that superfluid phonons are in general described by a scalar field 0 
governed by the effective field theory (EFT) [86] : 

C = P(X), X = 8-m<S> (4) 

where T is the gravitational potential. In particular, the type of superfluid, i.e., its equation of 
state, is uniquely encoded in the choice of P. 

Once we take seriously the idea that DM is a superfluid, the only question is — what kind of 
superfluid? The second central idea underlying this work is that DM phonons are described by the 
non-relativistic MOND scalar action, 


P(X) ~ AX^/\X\. 


(5) 


where A ~ meV to reproduce the MOND critical acceleration^] This choice corresponds to a 
particular superfluid, with P ~ p 3 . To mediate a MONDian force between ordinary matter, phonon 
must couple to the baryon density: 


Tint ~ 


A 

Mpi 


Opb ■ 


( 6 ) 


From a particle physics standpoint, such a coupling is fairly innocuous — it represents a soft explicit 
breaking of the global 17(1) symmetry. In the superfluid interpretation, however, where 9 is the 
phase of a wavefunction, this coupling picks out a preferred phase, which seems unphysical. One 
possibility is that ([6]) follows from baryons coupling to the vortex sector of the superfluid. This 
would give rise to a cos0pb operator [88l490| . thereby breaking the continuous shift symmetry down 
to a discrete subgroup. When expanded around the state at finite chemical potential 8 = pt, such 
operators would give ([6]) to leading order, albeit with an oscillatory prefactor. 

Thus, through ([5]) and (??), phonons play a key role by mediating a long-range force between 
ordinary matter particles. As a result, a test particle orbiting the galaxy is subject to two forces: 


3 The possible connection between MOND and superfluidity was mentioned briefly by Milgrom in [87] ■ We thank 
A. Kosowsky for pointing this out to us. 
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the (Newtonian) gravitational force and the phonon-mediated force. Our postulate is that the 
phonon-mediated force is MONDian, such that the DM superfluid reproduces the empirical success 
of MOND in galaxies. 

The fractional 3/2 power would be strange if ([5]) described a fundamental scalar field. As a 
theory of phonons, however, it is not uncommon to see fractional powers in cold atom systems. For 
instance, the Unitary Fermi Gas (UFG) [91, .92], which has generated much excitement recently in 
the cold atom community, describes a gas of cold fermionic atoms tuned such that their scattering 
length diverges USES]. The effective action for the UFG superfluid is uniquely fixed by 4d scale 
invariance at lowest-order in derivatives, £ufg(A") ~ A 5 / 2 , which is also non-analytic |95j |^| 

A hint on the nature of our condensate can be inferred from the (grand canonical) equation of 
state P{p), obtained by working at finite chemical potential 9 = fit: P ~ p 3 / 2 . Using standard 
thermodynamics, this implies a polytropic equation of state: 

P~p 3 . (7) 

We can compare this to the viral expansion P = k^T p + g- 2 {T)p 2 + gs{T)p 3 + ..., where the 
p term describes an ideal gas, the p 2 term describes 2-body interactions, the p 3 term 3-body 
interactions, etc. The P ~ p 3 dependence in our case suggests that DM particles have negligible 
2-body interactions and interact primarily through 3-body processes. It would be very interesting 
to find explicit examples of such superfluids in Nature and study in more detail their microphysical 
interactions. 

As is familiar from liquid helium, a superfluid at finite temperature (but below the critical 
temperature) is best described phenomenologically as a mixture of two fluids [98til00j : i ) the 
superfluid, which by definition has vanishing viscosity and carries no entropy; ii ) the “normal” 
component, comprised of massive particles, which is viscous and carries entropy. The fraction of 
particles in the condensate decreases with increasing temperature. Thus our framework naturally 
distinguishes between galaxies (where MOND is successful) and galaxy clusters (where MOND is 
not). Galaxy clusters have a higher velocity dispersion and correspondingly higher DM temperature. 
For m ~ eV we find that galaxies are almost entirely condensed, whereas galaxy clusters are either 
in a mixed phase or entirely in the normal phase. 

Assuming hydrostatic equilibrium with P ~ p 3 , the resulting DM halo density profile is cored, 
not surprisingly, and therefore avoids the cusp problem of CDM. Remarkably, for our parameter 
values (m ~ eV, A ~ meV) the size of the condensate halo is ~ 100 kpc for a galaxy of Milky-Way 
mass. In the inner region of galaxies where rotation curves are probed, the DM condensate has 
a negligible effect on baryonic particles, and their motion is dominated by the phonon-mediated 
MOND force. In the outer region probed by gravitational lensing, the DM condensate gives the 
dominant contribution to the force on a test particle. 

In the vicinity of individual stars the phonon effective theory breaks down and the correct 
description is in terms of normal DM particles. This is good news on two counts. First, it is well- 
known that the MONDian acceleration, while giving a small correction to Newtonian gravity in the 
solar system, is typically too large to conform to planetary orbital constraints. This usually requires 
introducing additional complications to the theory [101]. In our case, the MONDian behavior is 

4 Similarly, in the quasi-static limit ( 6 = 0) our action ~ X 3,72 becomes invariant under time-dependent spatial 
Weyl transformations: hij —> fl 2 (x,t)hij [96:1:97!. At lowest order in derivatives it is the unique action with this 
property. Intriguingly, the £0(4,1) global part of the 3d Weyl group coincides with the de Sitter isometry group, 
which hints at a deep connection between the MOND phenomenon and dark energy [97] . 
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avoided entirely in the solar system, as DM behaves as ordinary particles. The second piece of 
good news pertains to experimental searches of axion-like particles. By allowing the usual axion- 
like couplings to Standard Model operators, our DM particles can be detected through the suite of 
standard axion experiments, e.g., m- 

The superfluid interpretation has a number of observational consequences, discussed in detail in 
Secs. [9]—11H which can potentially distinguish this scenario from ordinary MOND and ACDM. We 
mention a few here: 


• As is well-known, a superfluid cannot rotate uniformly; when spun faster than a critical 
velocity, the superfluid instead develops localized vortices. The typical angular momentum 
of galactic haloes is well above the critical velocity, giving rise to an array of DM vortices 
permeating the galactic disc [TU31HTH] . Unfortunately these have negligible energy density, 
so their detection through gravitational lensing may prove challenging. Substructure lensing 
may soon be possible with the Atacama Large Millimeter Array [IQS] . 

• A key difference with ACDM is the merger rate of galaxies. Applying Landau’s criterion for 
superfluidity, we find two possible outcomes depending on the infall velocity. If the infall 
velocity is less than the phonon sound speed, then the galactic condensate halos will pass 
through each other with negligible dissipation. In this case the merger time scale will be 
much longer than in ACDM and involve multiple encounters, as dynamical friction between 
the superfluid halos will be negligible. If the infall velocity is greater than the sound speed, the 
encounter will drive halos out of equilibrium and excite DM particles out of the condensate. 
In this case dynamical friction will lead to a rapid halo merger, as in ACDM, and after some 
time the merged halo will thermalize and condense back to the superfluid ground state. 

• The story is even richer for merging galaxy clusters, such as the Bullet Cluster |106H108] . 
Here the outcome not only depends on the infall velocity, but also on the relative fraction 
of superfluid vs normal components in the clusters. If the infall velocity is sub-sonic, the 
superfluid components should once again pass through each other with negligible friction, 
however the normal components should be slowed down due to the significant self-interaction 
cross section. In general, we therefore expect that lensing maps of bullet-like systems should 
display two features: i) mass peaks coincident with the cluster galaxies, due to the (non¬ 
interacting) superfluid cores; ii) another mass peak coincident with the X-ray luminosity 
peak, due to the (interacting) normal components. Remarkably, this picture is consistent 
with the lensing map of the Abell 520 (MS0451+02) merging system [10914112] . The Bullet 
Cluster is also consistent with this picture if the sub-cluster (the “bullet”) is predominantly 
superfluid. 


The idea of a Bose-Einstein DM condensate (BEC) in galaxies has been studied before [103 . 1041 
I113)4l25j There are important differences with the present work. In BEC DM galactic dynamics 
are caused by the condensate density profile, similar to what happens in CDM, with phonons being 
irrelevant. In our case, phonons play a key role in generating flat rotation curves and explaining 
the BTFR. Moreover, the equation of state is different: the BEC DM is governed by two-body 
interactions and hence has P ~ p 2 , compared to ~ p 3 in our case. This difference only has a 
minor effect on the condensate density profiles, but it does imply a different phonon sound speed. 
In particular, for the Bullet Cluster the sound speed in BEC DM is only c s < 100 km/s, i.e., 


5 In the context of QCD axion, it has been argued that Bose-Einstein condensation can occur in galaxies 1126111271 . 
though this has been disputed recently 11281 . 
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more than an order of magnitude smaller than the bullet infall velocity. As a result dissipation is 
important, which puts BEC DM in tension with observations m- 


2 Dark Matter Condensation 


In order for DM particles to Bose-Einstein condense in galaxies, two conditions must be met. 
For the purpose of these initial estimates, we shall treat DM as weakly interacting particles for 
simplicity, leaving for future work a refined calculation including interactions. The first condition 

is that the de Broglie wavelength of DM particles AdB ~ ~ be larger than the mean interparticle 
/ \ 1/3 mV 

separation t ~ ( ^ J . This implies an upper bound on the mass: 


m < 



( 8 ) 


We shall apply this bound at virialization, which marks the initial moment when one can mean¬ 
ingfully talk about an individual halo. From standard collapse theory, virialization occurs when 
^ ~ 180. In terms of the present DM cosmological density — 3 x 10“ 30 g/cm 3 , the density at 
virialization is therefore 


p vir = (1 + z vir ) 3 180 p|^ ~ (1 + z vir ) 3 5.4 x 10 28 g/cm 3 . 


(9) 


Meanwhile, the velocity is related to the mass of the object as usual by H3Q] 

/ M \ 1 / 3 

Uvir = 127 1012^-1 jVf 0 J Vl + Zvir klll/s . (10) 

Substituting these into (|8j) , we obtain 

/ \ —1/4 

m < 2.3 (1 + ^ vir ) 3/8 y W 2 h -i M J eV - ( n ) 

Hence light objects form a BEC while heavy objects do not. Figure [2] shows the BEC region as a 
function of mass assuming z v ; r = 2 for concreteness. 

The second necessary condition for condensation is that DM particles thermalize, with the 
temperature set by the virial velocity. The interaction rate is given by [ 126} 


Afvpvir— , 
m 


( 12 ) 


where 


AC 


Pvir (27 t) 3 
m Ce (mu ) 3 


~ 10 3 (1 + 2 vir ) 3/2 


f m \ ~ 4 10 12 h l M P 


0 


M 


(13) 


is the Bose enhancement factor, i.e., of order 10 3 particles for a massive galaxy]^] The interaction 


rate should be compared to the dynamical time in galaxies, fdy 


Indeed, if the time scale 


a y n V&nPv„ 

for thermalization is shorter than the halo dynamical time, the coherence length for the condensate 
will be comparable to the size of the halo. This is necessary in order for phonons to act coherently 


Strictly speaking, (|12|) is valid provided that F <g mv 2 Il26i . which is easily satisfied in our case. 
















Figure 2: Dependence of the BEC region (shaded) on the DM mass m and the halo mass M. assuming 
z v ; r = 2 for concreteness. 


across the galaxy. Putting everything together, the condition Ttdyn ^ 1 can be expressed as a lower 
bound on the interaction cross section 


a 

m 


> (1 + z vir )~ 7/2 



M 


10 l2 h~ l M Q 


2/3 


52 


2 

cm 

g 


(14) 


Clearly the bound is most stringent for massive galaxies. Taking M ~ 10 12 h 1 Mq and assuming 
z v i r = 2 for concreteness, we obtain 


— > 
rsj 

m 



(15) 


We will see below that a mass of around 0.6 eV gives appropriate size halos, in which case 

__ 2 _ 

~ ^ 0.1 0 . • The lower end of this bound satisfies current constraints fT3THT33| on the cross 

section of self-interacting dark matter (SIDM) [134]. However, as we will see the phenomenology of 
superfluid DM is considerably different than SIDM, and each constraint much be carefully revisited. 

The resulting DM temperature is quite cold. The critical temperature can be readily obtained 
assuming equipartition, feeTc = \ m vh where v c saturates Q. The result is in the rnK range: 


T c = 6.5 



5/3 

(1 + z vir ) 2 mK . 


The temperature in a given halo, in units of T c , is 


T _ 0.1 / m \ 8 /3 / M 

¥ c ~ 1 + z vir VdV/ V10 12 h~ l M e 



(16) 


(17) 


At finite but sub-critical temperature, the system is phenomenologically described as a mixture of 
condensate and normal components. Neglecting interactions, the fraction of condensed particles 
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Figure 3: Fraction of DM particles in the condensate as a function of halo mass M for m = 0.4, 0.6 and 
0.8 eV, assuming ^ v ; r = 0. 


is [135] 


-Wond 

N 


= 1 - = 


T \ 


3/2 


~ 1 — 


0.03 


/ m \ 4 


M 


(1 + z vir ) 3 / 2 VeV/ 10 12 /i _1 Mq ’ 


T <T C . 


(18) 


Figure [3] plots the condensate fraction as a function of halo mass for z v j r = 0, for m = 0.4, 0.6 
and 0.8 eV. We see that galaxies (M < 10 12 h^ 1 M & ) are almost completely comprised of particles 
in the condensate, while massive clusters (lO 14 /i _1 M 0 <, M <, 10 15 /i _1 Mq) can have a significant 


fraction, if not all, of their particles in the normal phase. It is worth noting that (18) only holds 


for free particles; one expects the 3/2 power to change when including interactions. For instance, 
the power is 3 for particles trapped in a harmonic potential. We leave a careful calculation of the 
condensate fraction including interactions to future work. 


A few comments about cosmology. Since our DM particles are in the sub-eV mass range, they are 
axion-like particles. They must be produced out-of-equilibrium (e.g. through a phase transition) 
and remain decoupled from normal matter throughout the history of the universe. For instance, 
they can be generated through an axion-like vacuum displacement mechanism: in the early universe, 
the field is displaced from its minimum and starts oscillating once H < m. In this scenario DM 
particles are generated when Hi ~ m. The corresponding photon-baryon temperature is 


i; baryons ~ 


(19) 


which for m ~ eV is 50 TeV, i.e. around the weak scale! The velocity is initially relativistic, Vi < 1, 
and subsequently redshifts as v ~ 1/a. 

It is easy to see that, as soon as it generated cosmologically, DM becomes superfluid. Consider 
the de Broglie wavelength condition ([8]). Since v ~ 1/a and p ~ 1/a 3 cosmologically, both sides 
of the inequality are time-independent. Hence if Q is satisfied at any time it is satisfied at all 
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times. We can anchor this condition at matter-radiation equality using the observational constraint 
p eq ~ 10 -19 g/cm 3 ~ 0.4 eV 4 . Since v eq <C 1, it follows that 


m ~ 


P 1/A 
re q 


< 


fPeqV 74 

Wq/ 


( 20 ) 


hence the BEC condition is satisfied at all times. Similarly it is easy to show that thermalization 
proceeds efficiently, given the lower bound (15) on a jm and the high occupation number M 1. 


Naturally DM is much colder on cosmological scales than in collapsed structures. The temper¬ 
ature ratio T/T c = (v/v c ) 2 is constant cosmologically, where v c = (p/m 4 ) 1 / 3 saturates ([8]). Once 
again it is convenient to evaluate this at matter-radiation equality: 


2 /n^s/3 
eq VeV/ 


( 21 ) 


Assuming Vi ~ 1 when j\ bar y° ns ^ y/mMp\, we have u eq = , and therefore 


~ 10 


-28 ^/ W ^ 5 / 3 


( 22 ) 


which is very cold indeed. In contrast we see from © that T/T c ranges from 10 6 in dwarf galaxies 
(M ~ 10 6 Mq) to 10 -2 in massive galaxies (M ~ 10 12 Mq). In other words, cosmologically the 
DM superfluid can be described to an excellent approximation as a T = 0 superfluid. In collapsed 
structures, finite-temperature effects can be significant. As we will see, finite-temperature effects 
will be important in stabilizing the MOND phenomenon in galaxies. 


3 Superfluid Phase 


Once DM condenses and forms a superfluid, the relevant low-energy degrees of freedom are col¬ 
lective excitations in the form of phonons. Superfluid phonons are the Goldstone bosons for a 
spontaneously broken global U( 1) symmetry. In the non-relativistic regime they are in general 
described by a scalar field 9 with effective action [95] 


C = P(X ); X = 6 — mf - 


(Ml 

2m 


(23) 


where $ is the external gravitational potential, e.g., 4>(r) = £ or a spherical-symmetric 

static source. This effective Lagrangian is exact at lowest order in derivatives, with corrections 
suppressed by additional derivatives per field. To describe phonons at constant chemical potential p, 
we expand 

(Vc))) 2 


9 — pt -(- * 


X = p — + (j) ~ 


2m 


(24) 


In the case of interest, our conjecture is that the DM superfluid phonons are governed by the 
MOND action ([5]), 

P(I) = 2A(2 3 m)3/ ^ A-v^. (25) 

The square-root form is necessary to ensure that the action makes sense for time-like field profiles 
and that the Hamiltonian is bounded from below [69]. Note that the effective action (25) is only 
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well-defined away from X = 0, for both time-like and space-like profiles. In Sec. [6] we will give a 
more fundamental derivation of the phonon action starting from a complex scalar field with |<9'I'| 6 
interactions. As we will see, in that example a condensate only forms for 2m|X| > for some 
cutoff scale A r . 


To mediate a MOND force, phonons must couple to the baryon mass density pb- 

A 


Tint — ® 


Mpi 


6pb , 


(26) 


where ct is a dimensionless parameter. (The relativistic extension is more complicated and will be 
discussed in Sec. [8|) This operator explicitly breaks the shift symmetry only at the 1/Mpi level 
and is therefore technically natural. From the superfluid perspective, (26) can arise if baryonic 
matter couple to the vortex sector of the superfluid, giving rise to operators ~ cos Qpb that preserve 
a discrete subgroup of the continuous shift symmetry j88H90j . Expanding around a state at finite 


chemical potential, cj> = 6 — pt, this operator would yield a coupling of the form (26) with an 


oscillatory prefactor. For the purpose of the present work we shall treat (26) as an empirical term 
in our action necessary to obtain the MOND phenomenon. 

To summarize, our phonon theory depends on three parameters: the particle mass m. the scale 
A and the coupling constant a. The latter two parameters can depend on temperature, and thus on 
velocity, most naturally through the ratio T/T c . In particular they can assume different values on 
cosmological scales (where T/T c ~ 10~ 28 ) than in galaxies (where T/T c ~ 10~ 6 —10 -2 ). Specifically 
we will see in Sec. [7] that a must be ~ 10~ 4 smaller cosmologically, while A must be ~ 10 4 larger, 
in order to obtain an acceptable cosmology. The temperature dependence is therefore quite mild 
and can be ignored over the velocity range spanned by galaxies. Until Sec. [7] it will be implicitly 
understood that a and A assume their galactic values, ignoring any temperature dependence. For 
galaxy phenomenology, we will find in Sec. [4] that these two parameters must be related in order 
to reproduce the MOND critical acceleration: 


a 3 / 2 A = \J ao-Mpi ~ 0.8 meV 


a 


( A \- 2 / 3 

°' 86 ( meV ) 


(27) 


Hence a 0(1) for A rsj meV. 


3.1 Condensate and phonon properties 


The form of the phonon action (25) uniquely fixes the properties of the condensate through stan¬ 


dard thermodynamics arguments. We work at finite chemical potential, 6 = pt, setting phonon 
excitations and gravitational potential to zero. The pressure of the condensate is given as usual by 
the Lagrangian density, 

9 A 

(28) 


P(p) =—(2mpf/ 2 . 

This is the grand canonical equation of state P = P(p) for the condensate. Differentiating with 
respect to p yields the number density of condensed particles: 


n = = A(2m) 3 / 2 ^ 1 / 2 . 

Combining these expressions and using the non-relativistic relation p = mn, we find 

r? 

P = 


12A 2 m 6 


(29) 


(30) 
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This is a polytropic equation of state P ~ p 1+1 / n with index n = 1/2. In comparison, the standard 
DM BEC discussed in the literature is described by P ~ p 2 , corresponding to n = 1. We will see 
below that the halo profiles are nonetheless quite similar. 


Let us now consider phonon excitations on top of this condensate. Expanding (25) to quadratic 
order in phonon perturbations ft = 9 — pt, once again neglecting the gravitational potential, we 
obtain 


7-quad 


A (2m.) 3 / 2 
4/r 1 / 2 


ft 


2 — —— (V ft ) 2 

m, 


(31) 


The sound speed is 



(32) 


Expanding to higher order, we can identify the strong coupling scale of the theory. A typical 
interaction term is schematically of the form 


^higher-order D Am 3/ V /2 n d n </> r ' 


ft 


3/2 


1-1 


d n ft 


(33) 


where d stands for either dt or c s V, and the canonical variable is ft c ~ A 1 / 2 m 3 / 4 /r 1 ftft. The strong 
coupling scale, identified as the scale suppressing higher-dimensional operators, is 


A s = 



(34) 


3.2 Halo profile 


Given the equation of state (30), we can 
assuming hydrostatic static equilibrium, 
pressure and acceleration are related by 


compute the DM density profile of the condensate halo 
Focusing on a static, spherically-symmetric halo, the 


1 dP(r) 
p{r) dr 


d$(r) 

dr 

47tGn 

/2 



(35) 


Equivalently, since by definition p = mn = mj y this equation can be written as 


dX(r) d4>(r) 

dr = ~ m dr ’ 


(36) 


which automatically follows from the expression (24) for X when phonon excitations are set to 
zero. 

It is convenient to rewrite this equation in terms of dimensionless variables E and £, defined by 


p(r) 


P0 l 


jl/2. 


P0 


32-7rGNA 2 m 6 


e, 


(37) 
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Figure 4: Numerical solution to the n = 1/2 Lane-Emden equation with boundary condition s(0) = 1 
and S'(0) = 0. The solution vanishes at ~ 2.75. The dashed line is a simple approximate analytical form, 
3(£) =cos(f£). 


where po = p(0) is the central density. Differentiating (35) with respect to r, and expressing the 
result in the new variables, it is straightforward to obtain the n = 1/2 Lane-Emden equation^] 


J_ A ( _ _'=' 1/2 

e 2 d? v <) ~~ 


(39) 


The boundary conditions are 3(0) = 1 and 3'(0) = 0. The numerical solution, shown in Fig. [4j 
vanishes at 

6 - 2.75 , (40) 

which defines the size of the condensate: 


R = 


Po 


327rG , NA 2 m 6 


6 • 


(41) 


A simple analytical form that provides a good fit is 3(0 = cos (^^j, shown as the dashed 
curve in the Figure. The density profile is thus well approximated by 

p(r) ~ p 0 cos 1/2 (—) ; r < R. 

The central density is related to the mass of the halo condensate as follows [ 136 1 

M 6 


(42) 


Po = 


4vt/2 3 |H'(0)| ‘ 


' The Lane-Emden equation for general n is 


dS\ _ _„ n 
£ 2 d£ ^ d£j 

Analytical solutions exist for n = 0, 1 and 5 m- Other values of n require numerical integration. 


(43) 


( 38 ) 
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From the numerics we find H'(£i) ~ —0.5. Substituting (41), we can solve for the central density 

( M dm \ 2/5 (m\ 18 /s / A \ 6/5 _ 24 ( , 

(ev) (meVy 10 “ g/cm ' <44) 


Meanwhile the halo radius is 


R ~ 


/ M dm \ 1/5 /mx-6/5 / A \" 2 / 5 


V 10 12 Mp 


o 


VeV/ 


\meV J 


45 kpc. 


(45) 


Remarkably, for m ~ eV and A ~ meV we obtain DM halos of realistic size. In the standard 
CDM picture a halo of mass Mdm = 10 12 Mq has a virial radius of ~ 200 kpc. In our framework, 
the condensate radius can in principle be considerably smaller or larger depending on parameter 
values. For concreteness, in the remainder of the analysis we will choose as fiducial values 


m = 0.6 eV ; 


A = 0.2 meV . 


(46) 


(From (27) this corresponds to a ~ 5/2.) This implies a condensate radius of ~ 158 kpc for a halo 
of mass Mdm = 10 12 Mq. 

Through the relation p = mn = mjS, the above density profile fixes X(r ): 


X( r ) = 


P 


8A 2 m 5 
~ 2 x 


,„-e m dm V /5 / m V 1/5 f A ^ 2/! ’ ( 7Tr '\ ,n\ 

10 eV (lO0M0) (ev) (meV) C0S (2 R> ' (47) 


In particular, the central density determines the chemical potential: 


P = 


Po 


8A 2 m 5 ’ 


(48) 


which in turns determines the strong coupling scale (34): 


3/4 

A = _ p X _ 

S 8 3 / 8 A 1 /2 m 3/2 


~ meV ( J^) 3/W (^ ) 6/5 ( J_) 

- mcv ^ 10 12Mq) VeV/ \meVj 


2/5 


(49) 


Thus the strong coupling scale, like A, is of order meV. Finally, the gravitational potential <h(r) = 
mr 1 ( X(r ) — p) follows trivially from these relations. 

A few comments are in order. First we have neglected the effect of halo rotation in this cal¬ 
culation. Slowly-rotating BEC with polytopic equation of state can incorporated into a modified 


Lane-Emden equation j !37| . However we will see in Sec. 10 that rotating halos are typically unstable 
to the formation of quantum vortices, which carry the angular momentum. Second, R represents 
the size of the superfluid “core”, not of the entire halo. In reality we expect this core to be sur¬ 
rounded by DM particles in the normal phase, most likely described by a Navarro-Frenk-White 
(NFW) profile [138] . A careful analysis would require numerical simulations, which is beyond the 
scope of this paper. Third, the superfluid scenario offers a simple, if not mundane, resolution to 
the cusp-core and “too big to fail” problems H3 nsi- The density profile is cored and hence has 
a much lower central density than in collisionless CDM simulations, in better agreement with the 
inferred densities of MW dwarf satellites. 
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4 Including Baryons: Phonon-Mediated Force 


In this Section we derive the phonon profile in galaxies, modeling the baryons as a static, spherically- 
symmetric localized source for simplicity. We first focus on the zero-temperature analysis, where 
the Lagrangian is given by the sum of ( |25[ ) and ( |26| ). In this case we find two branches of solutions, 
depending on the sign of X. The branch with X > 0 has stable perturbations but does not admit 
a MONDian regime. The branch with X < 0 does admit a MONDian regime, where the phonon- 
mediated force approximates the MOND force law over the scales probed by galactic rotation curve 
observations, as desired. However, perturbations on this branch are unstable. Stability on the 


MOND branch can be restored by finite-temperature effects, as we will show in Sec. 4.2 


4.1 Zero-temperature analysis 

Recall our zero-temperature phonon Lagrangian: 


„ 2A(2m) 3 / 2 v r— A „ 

£=-3-xym-a—Op,,. 


(50) 


In the static spherically-symmetric approximation, 0 = pt + cj)(r), the equation of motion reduces to 

<VbH (51) 


V 


X\ V</> = 


2Afpi 

where X(r) = /i — m<h(r) — ^ . This can be readily integrated: 

^ . 


87rAfpp’ 2 


(52) 


The profile depends on the sign of X : 

• X > 0 branch: In this case the solution is 


cj)'{r) = x/m [p, - dp? - 2 


1/2 


m* 


where we have chosen the minus sign such that ft 
solution for X (r) is 


p = /! — m4> 
0 when 


-w-)=ha + ^-5 


(53) 

0. Equivalently, the 

(54) 


As a check note that X —> fi for Mb —> 0, which is consistent with our equation (36) for 


the density profile in the absence of baryons. More generally, we can solve (54) for the 

^2 


gravitational potential: [i = ft — m& = X + ^ v . Substituting into Poisson’s equation, we 
obtain 


V 2 [X + 


\ m 4 A f X \ 1/2 


4 m 2 X ) M 2 ! 


Pb 


m) + 2Mpi ’ 


(55) 


where we have used p = for the condensate matter density. In the absence of baryons, 


this reduces to the Lane-Emden equation (39). In the presence of baryons, it is easy to 


show that the solution is qualitatively similar, with the only notable difference being that the 
halo radius shrinks with increasing baryonic mass, as expected from the extra gravitational 
attraction due to baryons. 
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• X < 0 branch: On this branch the solution is 


1/2 


(f>'(r) = \fm ( fi + \j fi 2 + -^2 


m* 


(56) 


where we have dismissed a solution corresponding to imaginary cj)'. Equivalently, the solution 
for X(r) is 


1 


x ( r ) = o A - \/A 2 + —2 


ft 


777/ 


(57) 


Unlike the X > 0 solution, this branch admits a MONDian regime where k> /i, such that 

(58) 

In this limit the scalar acceleration on an ordinary matter particle is 


4,\r) ~ 


o *< r >= 12 


Mpi 


(59) 


To reproduce the MONDian result omond = \J we are therefore led to identify 

-2/3 


a 3 / 2 A = x/ooMpi ~ 0.8 meV ==> a — 0.86 ( - ] 

\ meV J 


(60) 


which fixes a in terms of A through the critical acceleration, as claimed earlier. That A is of 
order the dark energy scale is a direct consequence of the coincidence oq ~ Hq. 


Repeating the steps that led to (55), in this case we find 
VOX- " ' 


4 m 2 X) Mpj \ m ) 


+ 


Pb 

2Mpi 


(61) 


This equation generically leads to unphysical halos, with growing DM density as a function 
of r. The origin of this instability can be seen at the level of perturbations. Expanding (50) 
to quadratic order in phonon perturbations (p = cj> — 4>(r), we obtain 


Gquad — 


s,gn(X)A^(^-2^-2^( 

4y|X| V m m \ 


1/2 


\ 2X 


- - 2^- (X - z- - —Xdnp) 


2m J mr 2 


(62) 


The kinetic term ip 2 has the wrong sign for X < 0. 


To summarize, the X > 0 solution, given by (53), is continuously connected to the homogeneous 
condensate in the absence of baryons (Mb —> 0) and has stable perturbations. However, this 
branch does not admit a MONDian regime. The X < 0 solution, on the other hand, does admit 
an approximate MOND regime, but this branch has the peculiarity that <// remains non-zero even 
in the Mb —> 0 limit. Moreover, perturbations about this solution have wrong-sign kinetic term, 
indicating an instability. Below we will show that this instability can be cured by finite-temperature 
effects. 
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4.2 Finite-temperature effects 


The DM condensate in actual galactic halos has non-zero temperature, hence we expect that the 


zero-temperature Lagrangian (50) receives finite-temperature corrections in galaxies. At finite sub- 


critical temperature, the system is described phenomenologically by Landau’s two-fluid model: an 
admixture of a superfluid component, which has zero viscosity, and a normal component, which is 
viscous and carries entropy. The two components interact with each other. Their relative fraction 
is a function of temperature, and hence the mass of the collapsed object, as sketched in Fig. |3j 
At lowest order in derivatives, the effective field theory at finite temperature and finite chemical 
potential is [1391 

jC t ^ 0 = F(X,B,Y). (63) 


It is a function of three scalar quantities. The scalar X, already defined in (23), describes the phonon 
excitations. The remaining scalars are defined in terms of the three Lagrangian coordinates ^(x, t ), 
/ = 1, 2, 3 of the normal fluid0 


B = \J det d^ip 1 dVil> J ; 

Y = u 11 (p^O + mS ®) — m ~ fi — m4> + 0 + v ■ V (j >. 


(64) 


where = g tuKd^ 1 d^ K is the unit 4-velocity vector, and in the last step for Y 
we have taken the non-relativistic limit ~ (1 — 4>, v). By construction, these scalars respect the 
internal symmetries: i) ijj 1 —> ip 1 + cJ (translations); ii) ip 1 —> R 1 (rotations); in) ip 1 —> ^{ip ), 
with det = 1 (volume-preserving reparametrizations). 


Our goal is to seek a finite-temperature theory that will generate a MONDian phonon profile (58) 


over the scales probed by galactic rotation curve observations, while having stable perturbations and 
reasonable DM density profile. There is much freedom in specifying finite-temperature operators 
that will do the trick. The simplest possibility is to supplement (50) with the two-derivative 
operator 

(65) 


A£ = M 2 Y 2 = M' z ( jj, — m4> + 


where in the last step we have specialized to the normal fluid rest frame, v = 0. This leaves the 


static profile (56) unchanged, however it does modify the quadratic Lagrangian (62) by an amount 
A£q U ad = M 2 ip 2 . This will flip the sign of the kinetic term, and therefore cure the ghost, if 


M > 


Am 3 / 2 

vW 


0.5 


10 11 M & \ 1/4 


M b 


f—V /2 (—) 

\meV J \ 10 kpc J 


1/2 


m , 


( 66 ) 


which, remarkably, is of order eV ! Hence, for quite natural values of M, this two-derivative operator 
can restore stability. Furthermore, this operator gives a contribution A P = M 2 /r 2 to the condensate 


pressure, which obliterates the unwanted growth in the DM density profile mentioned below (61). 
Instead, the pressure is positive far from the baryons, resulting in localized, finite-mass halos. 

As another example, consider the finite-temperature Lagrangian: 


P(X,T) = 2AW!<1 l/ra _ M < 2 ’") 3/ A 


X — (3 (/x — + 


(67) 


8 In [139], Y is defined in terms of the relativistic phonon field 0 as Y = u^d^O. To translate to the non-relativistic 
description, we have substituted 0 = mt + 9 and subtracted the mass term. 
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where we have once again focused on the normal fluid rest frame. The dimensionless j3 parameter 
implicitly depends (mildly) on T/T c , though we will treat it henceforth as constant. This is of 
course a more ad hoc form of finite-temperature effects, but it has the advantage of facilitating the 
analysis. As we will see, in order to reproduce the MOND phenomenon with stable perturbations 
we will need 

( 68 ) 


'4 


in which case the quantity within absolute values is negative definite. 

First consider the DM density profile in the absence of baryons. Setting phonons and gravita¬ 
tional potential to zero, the pressure of the condensate is now given by 


P( / z,T) = ^0^(2m / r) 3 / 2 . 


(69) 


Thus the density profile is identical to the zero-temperature profile described in Sec. 3.2 modulo 
the replacement A — > ft/3 — 1A. For instance, instead of (45) the halo radius is now given by 

K<r) “ (i^%r (^r /5 (^v)" /5< ^- ir ' /54s kpc - <7o) 


Including baryons, the static, spherically-symmetric scalar equation becomes 


V 


+ 2m ( 3 1 j h -* \ _ ap h (r) 

Q> 1 “ 2 Mpi 


\Jft 2 + 2 m(f3 - 1 )/x 


(71) 


where ft = n — was introduced in (53). This integrates to 


k ' 2 + 2m ( ^ - 1 




ft ft 2 + 2 m(/3 - l)/x 


b' = n(r). 


(72) 


This implies a cubic equation for (j)' 2 , whose real root does not have a particularly illuminating 
analytic form. For concreteness we shall assume that /3 is strictly greater than 3/2. The solution 
then has the following behavior: sufficiently close to the baryon source, such that (j)' 2 ^$> mfi, the 
solution approximates the MOND profile (58), ft ~ ftn, and therefore scales as 1/r. Far from 


the baryons, such that ft 2 <C m/x, the solution tends to ft — \J which approximately 

scales as 1 jr 2 since jl is approximately constant. To summarize, assuming (3 > 3/2 the phonon 
profile is given by 

if r«r 4 , 

(73) 

if r>r*. 


k ~ 


3 

2 mp, 


(3-1 1 

TTn -o tv rv -' —tj 

2(3-3 r z 


The transition radius r* delineating these regimes occurs when k = mjl. It can be estimated by 
substituting the definition of k given in (|52j), and approximating fi as constant, with value set by 


the central density as in (48): /x ~ /9 q/ 8A 2 ?71 5 . The result is 


M b \ 1/10 /Mdm\" 2/5 /rnyS/5 / A \~ 8/15 


10 11 Mq 


V M b J VeV 


\meV ) 


20 kpc. 


(74) 
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Figure 5: Left Panel: The 0-mediated acceleration a$ (solid curve) is compared to the deep-MOND 
acceleration czmond (dashed curve) for a MW-like galaxy (Mb = 3 x 10 11 Af Q ) with cosmological DM-to- 
baryon ratio — 6, and fiducial values m = 0.6 eV and A = 0.2 meV. Right Panel: The ratio of 

the two accelerations as a function of radius is less than a few percent. 


For instance, with the fiducial parameters (46), the transition radius for a MW-like galaxy with 
Mb = 3 x 10 11 Mq and cosmic DM-baryon ratio ~ 6 is r* ~ 49 kpc. 

Figure [5] plots the numerical solution for <j>' , assuming (3 = 2 and the parameter values listed 
above. The Left Panel compares the scalar acceleration a<p (solid curve) to the MOND acceleration 
omond (dashed curve) as a function of r. The Right Panel shows the two accelerations only differing 
by a few percent, hence the predicted rotation curves are nearly identical to those of MOND. In 
particular, the “asymptotic” velocity is indistinguishable from that predicted by MOND (especially 
taking into account the uncertainties in the mass-to-light ratio), and the BTFR follows identically. 

It remains to compare the scalar acceleration a^ to the Newtonian acceleration odm due to the 
DM condensate profile. As we are about to show, in the MOND regime (r -C r*) the gravitational 
acceleration from the DM halo is negligible compared to the scalar-mediated MOND acceleration. 
In the opposite regime (r r*), on the other hand, the DM halo gives the dominant contribution 
to the force on a test baryonic particle. 

First, consider the regime r <C r* where a$ is approximately MONDian. In this case we have 
k 3> mfi, hence the DM density profile is 


A/2 


Pdm = {2mf /2 mAyJ/3 - 1 \/\X\ 
~ 2m 2 /3 — 1 \Jn(r ), 


(75) 


where we have made use of the substitution A —> y//3 — 1A mentioned earlier. Thus pdm ~ 1/r, 
and the Poisson equation can be straightforwardly integrated (ignoring baryons) to obtain acM = 
m Av/ ^~ 1 Y / /t(r)r 2 , which is constant. Comparing to the scalar acceleration a$ — y /k (r), we 


2Mp 

find 


«DM 

(Iff) 


yj (3 — 1 m 2 r 

2 a Mpi 


r 

0.6 — 

r* 


(r < r-*), 


(76) 


where in the last step we have assumed the parameter values listed below (74) for concreteness 


Hence, as claimed, the gravitational acceleration due to the DM halo is subdominant in the MON¬ 
Dian regime (r <C r*) and becomes comparable to the scalar-mediated acceleration around the 
transition radius r ~ r*. 
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Consider now the opposite regime, r r*. In this case we have X ~ p,, and the DM halo 
approximates the Lane-Emden density profile found in Sec. 3.2 The gravitational acceleration is 
(Idm = T \X'\ ~ while the scalar acceleration is a$ — , with q / ~ For the 

parameter values listed below (74) and taking (5 = 2 for concreteness, their ratio is given by 


QPM 

Cl(j) 


(r > r*). 


(77) 


Despite the crudeness of the estimate, this is remarkably consistent with (76) for r ~ r*. Hence, 
for r r*, the DM halo gives the dominant contribution to the acceleration on a test baryonic 
particle, as claimed earlier. 

Let us check the stability of the phonon background. Expanding (67) to quadratic order in 
perturbations tp = (f> — 4>(r), we obtain 


c 


quad 


A(2 m) 1 / 2 

Z 3 / 2 




(/3-1) 

(P- 1 ) 


(y-!)A + 

(f- 


2-1 

3 


</>' 2 ' 
2 m. 


1/ ■ / 

9W 


1 I A 2 + ~ ^ 


i/4 


2m 2 


— — i 
3 


- , ^' 2n 
M+ 2^ 


Z(dnipf 


9 


/ 2 


(78) 


- ^ 7/2 

where Z = (J5 — 1)A+^- The sign of the kinetic term is healthy if (5 > 1. Moreover, the sign 
of the (<9 q</?) 2 term is correct if (5 > 3/2, which ensures there are no gradient instabilities along 
the angular directions. Along the radial direction, the sign of the Lp' 2 is also correct if (5 > 3/2. 
It is then trivial to check by diagonalizing the kinetic matrix that radial perturbations propagate 
with the correct signature, he., they are free of ghosts or gradient instabilities. To summarize, the 
phonon background is perturbatively stable if (3 > 3/2, as claimed earlier. 

Note that, in the MOND regime (4>' 2 2 mfi), the phonon sound speed is c s ~ 4>'/m, which 
is enhanced compared to the sound speed ( [32] ) c s = \f2p/m computed in the absence of baryons. 
When we discuss various astrophysical probes below, we will nevertheless apply Landau’s criterion 
for the onset of dissipative effects using c s = y/2p/m, keeping in mind that this is conservative 
(since the actual sound speed is in fact larger). 


5 Validity of EFT and the Solar System 

_ i/2 

Our background solution (56) involves large phonon gradients, p, so naturally one should 

wonder whether it lies within the regime of validity of the EFT. First notice that in terms of the 
superfluid velocity v s = |V</>|/m and sound speed c s = \J2p/m, the MONDian regime ^ n 
precisely corresponds to v s 3> c s . It therefore violates Landau’s criterion v s < c s for the stability 
of superfluid flow. This is of course not surprising — Landau’s criterion is based on the stability 
of the superfluid against the creation of collective excitations, whereas we wish to work in a regime 
where baryons generate a large coherent phonon background. 
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As a check on whether this is legitimate, we can compare higher-derivative corrections to the 
leading-order P(X) Lagrangian. Such corrections by definition involve more than one derivative 
per field, hence they can be neglected as long as 


1 


A s 4>' A s r 


< 1 


(79) 


But since A s ~ meV, as we have seen in (49), this condition is trivially satisfied on astrophysical 
scales of interest. In other words, the phonon profile generated by a galaxy, while large relative to 
/r, is nevertheless very smooth on the scale of the cutoff. 

On the other hand, we should also verify that the local superfluid velocity does not exceed the 
BEC critical velocity, 

P \ 1/3 


v a < v c ~ 


m 


4 


(80) 


for otherwise the large phonon gradient will induce a local loss of coherence of the condensate. 
Equivalently, (80) can be understood as the requirement that the superfluid de Broglie wavelength 

i ( \ 1/3 

A ~ is much larger than the interparticle separation i ~ ^ J .To estimate v c , we use the 

halo mass density p = (2m) 3 / 2 mA-y/| X\ ~ 2m 2 A s /~K , where in the last step we have assumed the 
MOND regime. This gives 


0.025 


M b 


1/6 


10 11 M 0 J VeV 


(m\~ 2 / 3 / A \ 2//9 /kpc\ 


1/3 


\meV J 


\ r J 


(81) 


Meanwhile, the superfluid velocity is v s = cf>'/m ~ y/n/m, which gives 


v s ~ 0.008 


M b 


10 u M q J VeV 


1 ^ 2 /m\- 1 / A \ x ^ 3 kpc 


\ meV / 


(82) 


Thus the criterion (80) can be expressed as a bound on the distance from the galactic center: 

1/2 


r > 0.2 


M b 


lO n M 0 y VeV 7 


/ m \ - 1 / 2 / A \ 


-5/6 


meV J 


kpc. 


(83) 


This is satisfied down to the central regions of galaxies. 

The condition (80) does have important ramifications for the solar system. It is well-known 
within the standard MOND framework that the extra acceleration a albeit small compared to 
the Newtonian acceleration in the solar system, gives an unacceptably large correction to Newtonian 
gravity, in conflict with bounds from tests of gravity. One possible way out is to suitably modify 
P(X) at large X, but this requires fine-tuning [69]. Another possibility is to introduce a suitable 
higher-derivative galileon operator HDD, but this has the obvious disadvantage of complicating the 
theory. 

In our superfluid picture, we are naturally immune to this problem because the local phonon 
gradient generated by the Sun is so large that d80|) is violated throughout the solar system. Indeed, 


the superfluid velocity (82) due to the Sun (Afb = M 0 ) is 
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where r now represents the distance from the Sun. Meanwhile the BEC critical velocity (81) is 
set by the Milky Way galaxy (Mb = 3 x 10 11 Mq) evaluated at the location of the solar system 
(~ 8 kpc from the galactic center): 
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2/9 


(85) 


The criterion (80) can be expressed as a bound on the distance from the Sun: 
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( 86 ) 


which is larger than the solar system. 


The fact that (80) is violated in the solar system means that the BEC loses its coherence, and 


the condensate is replaced by a phase of normal DM. Hence the usual worries about MOND and 
local tests of gravity do not apply in our case. Furthermore, since our DM behaves as ordinary 
particles in the solar system, this is good news for direct detection experiments. By allowing the 
usual axion-like couplings to Standard Model operators, our DM particles can be detected through 
the suite of standard axion-like particle searches, e.g., m- 

6 A Relativistic Completion 

It is well-known that a superfluid can be described in the weak-coupling regime as a theory of a 
self-interacting complex scalar field with global 17(1) symmetry. The conserved charge associated 
with this symmetry is the total number of particles. A superfluid corresponds to a state which 
spontaneously breaks the global 17(1) and has finite charge density under this symmetry. 

In this Section we give an explicit example of such a theory that admits a condensate with 
P ~ /A 2 equation of state. After integrating out the radial mode, the resulting action for the phase 


to leading order in derivatives will be exactly given by (25), with the desired square root. The first 
theory that comes to mind is a scalar with hexic interactions, C = — |<1*1 2 — m 2 |<3?| 2 — A|<h| 6 . As 
shown in the Appendix, this gives P{X) ~ A' 3 / 2 , exactly the desired fractional power for MOND. 
However the sign is wrong. For a stable potential (A > 0), one is restricted to X > 0, hence spatial 
gradients can never dominate and the MOND regime is inaccessible. The MOND regime is only 
possible for A < 0, but this branch is of course unstable. 

Instead we will consider the following theory 
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(87) 


The scale A c is introduced to ensure that the theory admits a <h = 0 vacuum. The MOND 
regime corresponds to |<T>| 2 3> A 2 , as we will see shortly. Notice the absence of a quartic term 
((IA^>l 2 + m 2 |$| 2 ) . It is possible to include such a term provided its coefficient is not too large, 
as we will see towards the end of the Section. 

For our purposes it suffices to focus on the non-relativistic regime. Making the field redefinition 


T = pe ^+mt) ^ 


( 88 ) 
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and taking the non-relativistic limit, it is straightforward to arrive at 
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The power of p in the denominator of the second term guarantees the MOND scaling symmetries |96l 
EZI: assuming that spatial gradients dominate, and taking the MOND limit /)> A c , the action is 
invariant under the spatial scaling 
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The effective theory of the Goldstone mode is obtained by integrating out p. To leading order 
in the derivative expansion we can ignore (Vp) 2 contributions. In this limit the equation for p 
becomes algebraic: 


2m,X p 


= 0. 


_(A 2 + p 2 ) 7 + A 4 (2mX)V (A 2 — p 2 ) 

The MOND regime corresponds to p A c . Indeed, in this limit the solution is 

p 2 ~ Ay/2m (X 2 ) 1/4 = K^2m\X\ . 


Substituting this back into (89) gives, to leading order in derivatives, 

2A(2m) 3 / 2 


C ~ 


-x\/\x\- 


This agrees with the MOND phonon action (25). 
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The regulator scale A c implies that the MOND regime is restricted to p > A c , i.e., |X| > 2 c .$ . 


Using (59), this corresponds to 
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Observationally the MOND regime works quite well down to ~ ao/10, so this puts an upper bound 
on A c . By choosing A c a factor of a few smaller than A, the predicted breakdown could occur around 
the acceleration scale of the MW dwarf spheroidals, which are well-known to pose a challenge for 
MOND H3H17]. 

We can straightforwardly generalize the analysis to include a quartic term. To fast-track the 
discussion, let us immediately write the answer in terms of polar variables: 
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(95) 

where g is dimensionless. The power of p in the denominator of the new term is once again chosen 
such that (90) is a symmetry when A c and spatial gradients dominate. Focusing on this regime 

for simplicity, the equation of motion for p is a quadratic equation for p 4 . Choosing the branch 
such that the answer reduces to (92) as g —>• 0, we find 


p 2 ~ AW —gmX + 2m|X| ( 1 + 


2 \ V2 


(96) 


Upon substituting into (95), the action for the Goldstone will of course be different than (93), but 
will reduce to it in the limit of large X. What matters ultimately is that the Lagrangian for the 
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Goldstone has the same sign as (93), for both X positive and negative. It is easy to show that this 
is the case for 


5 <3- 


(97) 


Clearly the analysis can be generalized even further by including higher order terms, ((|<9 /i < f>| 2 + m 2 |<h| 2 ) n , 
with n > 4, provided they respect the scaling symmetry (90) in the appropriate limit. Their coef¬ 
ficients will be similarly constrained. 


7 Cosmology 


In this Section we study the cosmology of the DM superfluid. As mentioned towards the end of 
Sec. [2] the simplest genesis scenario is through a vacuum displacement mechanism, with DM being 
generated at a time when Hi ~ m corresponding to a baryon-photon temperature of order 50 TeV. 
The DM is initially very cold, it rapidly reaches thermal equilibrium with itself, but is decoupled 
from ordinary matter to first approximation. 


In order to obtain an acceptable background cosmology and linear perturbation growth, we will 
see that the A and a parameters of the phonon EFT must assume different values cosmologically 
than in galaxies. This is not unreasonable, as argued in Sec. [3], since these parameters are expected 
to depend on T /T c , and this ratio is ~ 22 orders of magnitude smaller cosmologically than in 
galaxies. Furthermore, we have already invoked finite-temperature effects in galaxies in Sec. 4.2 to 
ensure stability of the MOND regime. We will denote the cosmological values by Aq and a q . 


7.1 Equation of State 


The first thing to check is whether the condensate has sufficiently small pressure to act as dust. 
Recall from (30) our condensate equation of state: 


P _ p 2 

p 12Agm 6 


(98) 


The sound speed of linear fluctuations is identical, c 2 = w. At sufficiently low density (p <C Aom 3 ) 
the superfluid behaves as dust, whereas at high density (p 3> Aom 3 ) it behaves as a relativistic 
component. At the very least, we should impose that w < 1 at matter-radiation equality. Since 
w ~ 1/a 6 , and correspondingly c s ~ 1/a 3 , imposing w eq <C 1 will ensure that DM behaves to a 
very good approximation as dust throughout the matter-dominated era. Substituting the known 
value p eq ~ 0.4 eV 4 , this puts a lower bound on Aq: 


( m \ — 3 

—J eV. (99) 

In particular Ao 3> 0.5 eV for our fiducial value m = 0.6 eV. This is roughly four orders of magnitude 
larger than the fiducial value A = 0.2 eV assumed in galaxies. This can be achieved, for instance, 
if A depends on temperature as 
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7.2 Coupling to baryons 


The above equation of state was derived ignoring the coupling to baryons. We now rectify this and 
derive the phonon cosmological evolution sourced by the baryonic density. Setting 6 = 9(t ), the 


phonon action given by (25) and (26) becomes 
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2Aq (2m) 3 / 2 a 3g3/2 
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Varying with respect to 6 gives the equation of motion 
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Since a 3 pb = const., we can integrate straightforwardly: (2m) 3 / 2 # 1 / 2 = —where C 
is an integration constant. In the non-relativistic approximation, the energy density is p = mn = 
mAo(2m) 3 / 2 0 1//2 , hence 

P = ~°YY^ mt Ph + Pdust, (103) 

Mpi 

where Pdust = mApC /a 3 . This term is recognized as the dust contribution studied in the (baryon- 
free) analysis of Sec. 7.1 Note that the non-relativistic approximation breaks down when 6 ~ m, 


corresponding to p ~ m 3 Ao, which from (98) is precisely when pressure becomes non-negligible. 

In the matter-dominated era, t ~ a 3 / 2 , the baryonic contribution ~ pbf redshifts as 1/a 3 / 2 
whereas the second term redshifts as usual as Pdust ~ 1/a 3 . In order for the superfluid to behave 
as ordinary dust, the second term should dominate over the first all the way to the present time: 


apAo , Pb 
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(104) 


Substituting the age of the universe to = 13.9 x 10 9 yrs ~6x 10' 32 eV 1 , and assuming a DM-to- 
baryon ratio of Pdust/Pb = 6, we obtain 
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(105) 


where the last step follows from (99). In particular, ao 10 4 for our fiducial value m = 0.6 eV. 


This is roughly four orders of magnitude smaller than the value a = 2.5 obtained in galaxies by 
matching to MOND. This can be achieved, for instance, if a depends on temperature as 


a(T) = a 0 (l + K a (T/T c 


K n ~ 10 4 . 


(106) 


Note that while A(T) and ct(T) both depend on temperature, the scale A' ~ aA appearing in 
the phonon-baryon coupling (|26l) is nearly temperature-independent. 


7.3 Velocity-dependent critical acceleration 

An immediate corollary of A(T) and a(T) being temperature-dependent is that the critical accel¬ 
eration, 
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also depends on temperature. More precisely, since the product aA is constant to a first approx¬ 
imation, the temperature dependence of ao is governed by a. In particular, in light of (105) we 
obtain 

4 ao, (108) 


a c 0 osmo < 10 


- 4 , 


where ao is the typical MOND value ([2]) in galaxies. Given this strong suppression of ao, it follows 
that gravity is highly Newtonian on cosmological scales. 

Another consequence is that there is no longer a universal value for the MOND critical accel¬ 
eration in galaxies, instead ao is predicted to depend on the velocity dispersion. The functional 
dependence is model-dependent of course, but the generic trend is that ao decreases with decreasing 
velocity. Intriguingly, this trend has been noted in the data — low-surface brightness galaxies tend 
to prefer a lower value of ao m- 


8 Gravitational Lensing 


In the context of TeVeS [59], the absence of DM in galaxies forces one to assume a rather complicated 
coupling between the scalar field 4> and matter fields in order to reproduce acceptable gravitational 
lensing. For starters, one supplements the theory with a 4-vector field A fl , which is unit time¬ 
like g^ v A^A U = —1. Then the non-relativistic scalar-matter interaction £ CO upimg = 
is covariantized by coupling matter fields to an effective metric gj^ , defined in terms of the 
Einstein-frame metric g^ via 
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(109) 


In the weak-field, quasi-static regime, g^ takes the usual form: goo = — (1 + 2$), go* = 0 and 
gij = (1 — 2 &)6ij. To this order we can ignore perturbations in the vector field, i.e., A M = (1, 0, 0,0), 
such that 
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where is of course sourced by baryons only: 
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This line element is exactly of the General Relativity form, albeit in terms of a shifted gravitational 


potential T + Hence the mass inferred from lensing observations matches the mass inferred 


A/pi 


from dynamical measurements. The TeVeS metric (109) was of course precisely engineered for this 
purpose. Specifically, the equality of gravitational potentials in (110) traces back to the precise 


factor of 2 in the combination g tlu + 2A^ l A v appearing in (109). This relative factor is not protected 
by any symmetry. 

In our case the story is simpler on two counts. First, there is no need to postulate an additional 
vector field. The normal fluid component already provides us with a time-like vector field u^, as 


discussed in Sec. |4.2| Second, the DM in galaxies contributes to lensing, hence the TeVeS factor of 
2 can be generalized: 

9ii,2 - 9fiu ~ 9\iv + (1 + 7 , ( 112 ) 
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with 7=1 corresponding to the TeVeS tuning. Working in the rest frame of the normal fluid, this 
gives in the weak-field limit, 
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where $ is now sourced by both baryonic and dark matter: 

V 2 <f> = 47 tGn (pb + Pdm) • 


(113) 


(114) 


Hence the lensing signal will arise from a combination of the 7 term in (113) and the DM con¬ 


densate density profile shown in Fig. [4} Determining the allowed range of 7 will require a detailed 
comparison with lensing observations, which is beyond the scope of this paper. What is clear is 
that there should be considerably more freedom than in TeVeS. It may even be that 7 = —1 is 
allowed, in which case the coupling to matter would reduce to a simple conformal coupling. 


In most of our discussion so far, we have assumed fiducial parameter values (46) such that the 
condensate radius is of order the virial radius, e.g. R 158 kpc for Mdm = 10 12 M 0 compared 
to 200 kpc for the virial radius. By choosing other parameter values, however, we can consider 
smaller condensate radii, in which case the condensate core will be surrounded by an envelope of 
DM particles in the normal phase, presumably with a NFW density profile. In that case the lensing 
signal could result primarily from the NFW envelope. This deserves a dedicated analysis, which 
will appear elsewhere. 


9 Merging Clusters: the Bullet and the Counter-Bullet 


The “Bullet” Cluster 1E0657-57 [106H108] shows lensing peaks displaced from the gas and centered 
around the galaxy distribution. This is expected in CDM: the halos are made up of weakly inter¬ 
acting dark matter particles that fly past each other, together with the galaxies, while the baryonic 
plasma is slowed down by ram pressure and ends up spatially segregated from the halos. By now 
observers have identified over thirty such merging systems [851 . 1141] . 

Galaxy clusters in the present context are comprised, either partially or fully, of DM particles 
in the normal phase. Hence we also expect lensing peaks displaced from the gas, due to the DM 
component. An important consideration is the constraint this imposes on the self-interaction cross 
section of the DM [133 . [142 ]. The tightest constraint comes a recent analysis of ~ 30 merging 
systems 
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(115) 


At face value there is a window for which this is consistent with our lower bound (15) for DM 


condensation in galaxies. However we think that the constraint (115) is not as stringent in our case. 


Indeed, (115) was derived assuming a single DM component, whereas the 2-fluid mixture makes for 


a much richer situation. The heterogeneous nature of merging systems, with different interactions 
among their components, can result in a significantly weaker bound in our casej^] Specifically, 
we expect the superfluid components to pass through each other with negligible dissipation if the 
relative velocity is sub-sonic, 

+nfall Cs ■ ( 116 ) 


9 This loophole was also exploited recently with ultra-strongly interacting DM ns]. 
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Using (32) and (47), and assuming the fiducial parameter values (46) for concretenes, it is straight¬ 
forward to show that c s — 1400 krn/s for the sub-cluster (M su b — 10 14 Mq), while c s — 3500 krn/s 
for the main cluster (M ma i n ~ 1O 15 M0), assuming a significant fraction of their mass is condensed. 


These values are comparable to the estimate of ~ 2700 krn/s for the relative velocity (144( 114'5] , 
indicating that dissipative processes between the superfhhd cores should be suppressed!* 3 ] 


In general our framework predicts two distinct features that should appear simultaneously in the 
lensing maps of bullet-like merging systems: i ) mass peaks coincident with the cluster galaxies, due 
to the (non-interacting) superfluid cores; ii ) another mass peak, approximately coincident with the 
X-ray luminosity peak, due to the (interacting) normal components. Interestingly, this is consistent 
with the complex mass structure of the “train wreck” Abell 520 (MS0451+02) merging system (109b 
ra, often hailed as a counterexample to the “Bullet” cluster. Aside from the “bullet-like” lensing 
peaks around bright galaxies segregated from the gas, this system also exhibits a puzzling “dark 
core” overlapping the X-ray gas without corresponding bright galaxies. In the context of SIDM, 
the cross section required to explain this feature is inconsistent with the bullet bound ( |115[ ) |112j . 
In our case, however, the dark core is naturally explained as due the normal DM components 
Intriguingly, even in the case of the Bullet Cluster the combined strong and weak lensing map 
reveals a significant mass peak coincident with the X-ray gas [108j . 


Another way that (15) and (115) can be satisfied simultaneously is if the cross section is velocity- 
dependent. This is in fact expected for dark atoms, since the cross section between ordinary 
atoms is generally a rich function of velocity [ 147j due to various atomic resonances. Such velocity 
dependence may imply a suppressed cross section in clusters, where the typical virial velocity is 
~ 10 times larger than in galaxies. A velocity-dependent cross section was proposed in the SIDM 
context to simulateously match the inferred profiles of dwarf galaxies and galaxy clusters I148H152] . 


10 Vortices 

As is well-known, a superfluid cannot rotate uniformly. When spun faster than a critical angular 
velocity the superfluid develops quantum vortices that carry the angular momentum [ 1531 . In the 
context of BEC dark matter, vortex formation was initially considered in [ T03j and studied in 
detail subsequently in (154] . For the purpose of this paper, we shall content ourselves with simple 
dimensional analysis along the lines of [103] , 

We can immediately convince ourselves that our halos rotate much faster than critical velocity. 
The critical angular velocity for vortex formation in a vessel of radius R is, up to a logarithm 
factor (153j . 

“ cr ~ ~ 10 ' Jls "‘ ■ (117) 

where we have assumed a halo radius R ~ 100 kpc and mass m ~ eV. On the other hand, 
the angular frequency of a DM halo of density p is u> ~ Xy/G^p, where A = J^ M 5/2 is the so- 
called spin parameter, while L and E are the total angular momentum and energy of the halo 
respectively. From CDM simulations one finds 0.01 < A < 0.1. Substituting a typical density of 
order p ~ 10~ 25 g/cm 3 , we find 

uj ~ 10 _18 A s -1 . (118) 

10 This is unlike BEC DM, where the critical velocity is only ~ 100 km/s 11291 

11 It has been argued that the contradictory nature of the Bullet and counter-Bullet can also be explained in the 
BEC DM context |146j . 
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Hence u S> w CT ; and vortex formation is unavoidable. 

The line density of vortices can also be readily estimated, 


<t v ~ moj ~ 10 2 A AU 2 . (119) 

In a galactic halo of radius R ~ 100 kpc, this means N v ~ 10 23 vortices in total. Their core radius 
is of order the healing length £, which is estimated as 


£ ~ -~ mm , 

mc s 


( 120 ) 


where we have assumed a halo of mass M ~ 10 12 Mq and used the fiducial parameters (46). Thus 
the core radius is an order of magnitude or so larger than the average interparticle separation in 
galaxies. 


It would be interesting to study whether these vortices can be detected observationally, for 
instance through gravitational lensing. This may prove challenging, since their kinetic energy per 
unit volume is tiny: A p ~ —p ~ 10 _ 33 Ap. Substructure lensing may soon be possible with the 
Atacama Large Millimeter Array ra¬ 


il Other Astrophysical Consequences 

In this Section we speculate on various astrophysical implications of superfluid DM. For the purpose 
of this initial paper our discussion will be quite qualitative, leaving a more careful analysis to the 
future. 


Galaxy mergers: A very interesting question is what happens during galaxy mergers. Fol¬ 
lowing Landau’s criterion for superfluidity, the merger dynamics depend on the infall velocity 
Uin fa.i i compared to the phonon sound speed c s within halos. The sound speed in a given halo 
is generally of order of the virial velocity. For instance, for our fiducial parameter values (46) 
we find c s ~ 220 krn/s in a 1O 12 M 0 halo. If the infall velocity is ultra-sonic, Uinfaii > c s , the 
encounter will drive halos out of equilibrium, exciting DM particles out of the condensate. 
As in ACDM, dynamical friction will lead to a rapid halo merger, and after some time the 
merged halo will thermalize and condense back to the superfluid state. If the infall velocity 
is sub-sonic Ui n f a u < c s , on the other hand, the merger time scale will be much longer and 
involve multiple encounters, as dynamical friction between the superfluid halos will be negli¬ 
gible. This is similar to what happens in MOND 


• Reduced dynamical friction: The overall reduction in dynamical friction due to the su¬ 
perfluid nature of the DM halo alleviates a number of minor problems with CDM. Instead 
of being slowed down by dynamical friction, galactic bars in spiral galaxies should achieve a 
nearly constant velocity, as favored by observations ra- This effect has been pointed out in 
BEC DM [1171 1156] and MOND [42j . Reduced dynamical friction would also help with the 
M81 group of galaxies — see and references therein. 

Another interesting system is the Fornax dwarf spheroidal]^] Five satellite globular clusters 
orbit Fornax close enough that they should lie within their host’s DM halo, assuming an 

12 We thank Lam Hui for pointing this out to us. 
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NFW profile. If so, however, dynamical friction should have caused the globular clusters to 
rapidly fall towards the center of Fornax [ 1581 l"59j ] . In reality Fornax shows no sign of such 
mergers. A possible explanation in ACDM is that the Fornax’s DM halo is cored, with the 
globular clusters orbiting on the periphery [160]. In our case, the situation is unclear, due 
to two competing effects. On the one hand, dynamical friction within Fornax’s superfluid 
DM halo should be reduced, as already mentioned. On the other hand, dynamical friction 
with stars is enhanced in MOND, thereby reducing the merger time USB- This will require 
a detailed study. 

• Dark-bright solitons: Given the large coherence length of the BEC, galaxies in the process 
of merging should exhibit interference patterns (so-called dark-bright solitons) that have been 
observed in counterflowing BECs at super-critical velocities, e.g., m- This effect has been 
studied to some extent in ultra-light BEC DM [163]. It would be interesting to estimate the 
spatial extent and lifetime of the fringes to see whether they are potentially observable. It is 
intriguing to speculate that this can offer an alternative mechanism to generate the spectac¬ 
ular shells seen around elliptical galaxies [164] 


• Vast planar structures and tidal dwarfs: The vast planar structures seen in the Local 
Group P3H23] and beyond [2lJ find a possible explanation in our scenario, similar to that 
proposed in MOND [TSJ - Namely, the planar structures around the MW and Andromeda 
would be the result of tidal stripping during a fly-by encounter between these galaxies. In 
particular, most of their satellite galaxies would be tidal dwarfs. With the MOND force law 
it has been estimated that MW and Andromeda had a fly-by encounter ~ 10 Gyr ago, 
with < 55 kpc closest approach distance [39]. In ACDM, such a past encounter, while 
in principle possible, would have disastrous consequences: dynamical friction between the 
extended halos would cause a rapid merger of MW and M31. In MOND, however, there is 
only stellar dynamical friction and a merger can be avoided ]40l442] . 

Similarly, in our case dynamical friction is suppressed among DM particles if the infall ve¬ 
locity is sub-sonic, as mentioned before. If even a tiny amount of superfluid DM is stripped 
along with the tidal dwarf galaxies created in the process, their dynamics will be governed by 
MOND, resulting in flat rotation curves that fall on the BTFR, consistent with observations 
of the NGC5291 dwarfs [ 321 [33] . 

• Globular clusters: It is well-known that globular clusters contain negligible amount of DM. 
Indeed, their observations are well-fitted by taking only the baryonic mass into account and 
assuming Newtonian gravity. This poses a problem for MOND |52j . Our case is clearly 
different, since the presence of a significant DM component is necessary for the MOND phe¬ 
nomenon to occur. To the extent that whatever mechanism (e.g., tidal stripping) responsible 
for DM removal in ACDM is also effective in our case, our model predicts DM-free (and 
therefore MOND-free) globular cluster dynamics. 

• Tri-axial DM halos: A key prediction of collisionless CDM simulations is the ellipticity 
of DM halos [165 ], which is borne out by lensing observations. Lensing mass reconstruction 

1J We thank Ravi Sheth for suggesting this idea to us. 
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of galaxy clusters often require an elliptical DM clump around the brightest central galaxy. 
On the other hand, DM self-interactions tend to isotropize the DM distribution, resulting 
in more spherical halos. To match the ellipticity of galaxy cluster MS2137-23 inferred from 
strong lensing observations, m claimed an even tighter bound than ( |115[ ), though recent 
SIDM simulations find consistent halo morphology for cross sections as large as ~ cm 2 /g |166j . 


Since superfluids have surface tension, the superfluid core surrounding the brightest central 
galaxy should be highly isotropic. The source of ellipticity must be the subdominant normal 
DM component. The normal-normal self-interaction cross section ~ 0.1 cm 2 /g is consistent 
with the observational bound [ 166] . However, since the normal component only makes up 
a small fraction of the total DM mass in the central region of galaxy clusters, the rate of 
self-interaction is considerably smaller, and much larger cross sections are therefore allowed. 
This clearly deserves further study. Interestingly, the ellipticity has been observed to decrease 
towards the center of clusters (r < 16 kpc) [167 ]. consistent with a highly spherical superfluid 
core. 


12 Discussion 

In this paper we proposed a novel theory of DM superfluidity that reconciles the stunning success 
of MOND on galactic scales with the triumph of the ACDM model on cosmological scales. The DM 
component consists of self-interacting axion-like particles which are generated out-of-equilibrium 
and remain decoupled from baryons throughout the history of the universe. Provided that its mass 
is sufficiently light and its self-interactions sufficiently strong, the DM can thermalize and form a 
superfluid in galaxies, with critical temperature of order ~mK. The superfluid phonon excitations 
are assumed to be described by a MOND-like action and mediate a MONDian acceleration on 
baryonic matter. Superfluidity only occurs at sufficiently low temperature, or equivalently within 
sufficiently low-mass objects. This naturally distinguishes between galaxies (where MOND is suc¬ 
cessful) and galaxy clusters (where MOND is not): due to the larger velocity dispersion in clusters, 
DM has a higher temperature and hence is either in a mixture of superfluid and normal phase, or 
fully in the normal phase. 

The superfluid interpretation makes the well-known non-analytic nature of the MOND scalar 
action much more natural. The phonons of the unitary Fermi gas, which has attracted much 
excitement in the cold atom community recently m, are also governed by a non-analytic kinetic 
term (with 5/2 power instead of 3/2 for our DM superfluid). The DM condensate equation of state 
P ~ p 3 suggests that our superfluid arises from three-body interactions. It would be fascinating to 
find precise cold atom systems with the same equation of state as our DM condensate. Practically 
this would yield important insights on the microphysical interactions that give rise to this particular 
superfluid. Tantalizingly, it might allow laboratory simulations of the properties and dynamics of 
galaxies. 

The rich physics of superfluidity leads to a number of observational signatures that can poten¬ 
tially distinguish our scenario from ordinary MOND and/or standard ACDM: numerous low-density 
vortices in galaxies; merger dynamics depending on the infall velocity vs phonon sound speed; dis¬ 
tinct mass peaks in bullet-like cluster mergers, corresponding to superfluid and normal components; 
interference patters in super-critical mergers. Studying these observables with numerical simula¬ 
tions promises to be fascinating. 
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Appendix: Why $ 6 Fails to Give MOND 


In this Appendix we show that a complex scalar held with hexic interactions yields a phonon action 
with the desired 3/2 power but with the wrong sign to give the MOND phenomenon. Our starting 
point is the relativistic action 


C = -|d„T| 2 -m 2 |T| 2 - -|T| b , 


(A-I) 


This theory is invariant under global 17(1) symmetry, with associated conserved charge being the 
number of particles. Making the replacement T = &e imt and taking the non-relativistic limit, the 
theory becomes 


C = - (T<9 t T* - T^T) - - —\|T| 

2 v ’ 2m 24m 31 1 

The equation of motion is a nonlinear Schrodinger’s equation, 


- idt'b + 


V 2 T 


m 


8m 3 


]T| 4 T = 0. 


(A-II) 


(A-III) 


This equation possesses the following homogeneous background solution which describes the BEC 
at zero temperature 

T 0 = V2mve l,lt , (A-IV) 

\ 4 . 

where p = the chemical potential. Meanwhile v is related to the number density of particles 

in the condensate, n = 2 mv 2 , which in turn is the Noether charge density of the spontaneously 
broken 17(1) symmetry. 

To study the spectrum of perturbations around (A-IV), we can expand as follows 

T = y/2m(v + p)e^ t+(p) , (A-V) 

where p is the perturbation of the order parameter, while </> is the Goldstone bosonj^] Substituting 
into (A-II) we obtain 

C = —(V/?) 2 + 2 m{y + p) 2 


p + cf)- 


(V^ 

2 m 


-7(» + rt«. 


(A- VI) 


4 Strictly speaking, (A-IV I spontaneously breaks the diagonal combination of the internal (7(1) and time translation. 
Therefore (j) is the Goldstone boson of this symmetry. 
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The low energy spectrum of the theory can be deduced as usual by linearizing the equations of 
motion and computing the characteristic determinant. The analysis shows there is one dynamical 
degree of freedom in the spectrum, with dispersion relatiorp*] 

w 2 = k 2 + 0(k 4 ). (A-VII) 

rn z 

Thus A > 0 is necessary for stability. 

The effective theory of the Goldstone can be obtained by integrating out p. To leading order in 
the derivative expansion the (Vp) 2 term can be ignored, with resulting action 


£ 






(WO 2 

2 m 


1/2 


(A-VIII) 


As a consistency check let us linearize the theory and compare the result to the dispersion rela¬ 


tion (A-VII). The quadratic Lagrangian for reduces to 


_ (2mp\ l/2 / 1 - 2 1 . - ,, 2 

^quad — 777-1 ^ I ( 0 0 ^ v^ 7 


2 (1 


777 - 


(A-IX) 


Perturbations are stable for p > 0, which is guaranteed by A > 0. Taking into account the explicit 


expression (A-IV) for the chemical potential we recover the dispersion relation (A-VII) 


Notice that (A-VIII) looks very similar to (25). It involves the correct fractional power needed 


for the MOND action. Unfortunately, because of the requirement A > 0 the gradient term can 
never dominate over p, and the would-be MOND regime is inaccessible. One may be tempted to 
focus on A < 0 instead, since in that case the limit of large gradients appears to be well defined. 
Moreover, we even obtain the correct equation of state for the condensate when we set = 0, 


taking into account that p/X > 0. However, according to (A-IX) the perturbations around the 
condensate have a ghost-like kinetic term for p < 0. The physical origin for this instability is very 
simple — A < 0 corresponds to an attractive interaction between bosons, hence the homogeneous 
BEC is unstable against collapse. 

In contrast the theory with |9<I>| 6 interactions studied in Sec. [6] precisely gives the phonon 


theory (25) and has stable perturbations around the homogeneous BEC background. 
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